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ABSTRACT 


We  report  on  R&D  that  is  enabling  enhancements  to  the  Bayesian  Infrasound  Source  Location  (BISL)  method  for 
infrasound  event  location.  The  focus  of  this  effort  is  on  improving  BISL  through  the  development  and 
implementation  of  physics-based  priors.  Phase  identification  in  BISL  is  incorporated  through  a  prior  probability 
density  function  on  group  velocity.  We  report  on  the  development  of  a  distance-dependent  prior  through  a  massive 
simulation  effort  that  will  be  validated  by  a  ground-truth  data  analysis  effort.  In  addition  to  providing  constraints  for 
the  prior  Probability  Density  Function  (PDF),  this  coupled  effort  is  also  providing  constraints  on  the  probability 
distributions  that  are  used  in  the  likelihood  equations  in  BISL.  We  assess  the  enhancement  in  location  precision 
using  a  dynamic  prior,  instead  of  the  previous  uniform  prior,  for  some  example  events.  Finally,  we  discuss  efforts  to 
improve  the  numerical  implementation  of  BISL  using  Monte  Carlo  integration. 
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OBJECTIVES 


This  research  aims  to  enhance  infrasound  event  localization  by  utilizing  both  empirical  observations  and  model- 
based  propagation  constraints  in  the  BISL  algorithm  (Modrak  et  ah,  2010)  in  the  form  of  prior  information. 

The  specific  objectives  of  this  research  are  discussed  below: 

Modeling  infrasound  propagation  under  a  comprehensive  set  of  atmospheric  scenarios. 

To  first  order,  vertical  variations  in  atmospheric  temperatures,  wind  speeds,  and  wind  directions  determine  the 
effective  velocity  and  thus  propagation  paths  at  regional  and  global  distances.  To  capture  the  variability  of  acoustic 
propagation  in  time  and  space,  multi-year  multi-scale  atmospheric  specification  datasets  are  used  to  run  the 
numerical  modeling  with  rays  propagating  at  representative  directions  and  elevation  angles  covering  regional  and 
global  distances. 

Extracting  information  from  modeling,  and  generating  and  implementing  prior  information  into  BISL. 

By  modeling  infrasound  propagation,  we  can  classify  ray  paths  that  return  to  ground  with  the  respective  velocities, 
ranges,  and  maximum  elevation.  This  information  leads  to  the  generation  of  propagation  catalogs  that  can  be  used 
for  the  generation  of  prior  information  that  is  fed  into  the  localization  algorithm  to  constrain  the  solution  to 
physically  meaningful  scenarios. 

Estimating  accuracy  and  precision  for  standard  and  enhanced  prior  PDF  for  BISL  implementation. 

Measuring  the  successfulness  of  the  enhanced  localization  requires  the  introduction  of  metrics  that  account  for  the 
precision  and  accuracy  of  the  solution.  The  BISL  localization  algorithm  is  applied  to  global  and  regional  ground- 
truth  events  to  estimate  the  95%  credibility  area  and  the  offset  between  estimated  and  real  events  to  assess  the 
precision  and  the  accuracy  of  the  algorithms. 

Exploring  further  enhancements  for  localization  schemes. 

Our  planned  mathematical  manipulation  of  the  BISL  framework  will  allow  further  enhancements  by  allowing  the 
likelihood  density  function  to  probe  more  complex  parameter  spaces. 

RESEARCH  ACCOMPLISHED 


Overview  of  the  BISL  and  Theory  of  Enhancement. 

The  BISL  method  estimates  the  most  likely  source  location  and  origin  time  by  performing  a  grid  search  over  a  4D 
parameter  space  (Modrak  et  al.,  2010).  This  multidimensional  space  includes  the  two  Cartesian  coordinates  ( xQ  and 

jpQ),  velocity  (V),  and  time  ( t0).  In  the  Bayesian  framework,  the  posterior  probability  density  function  (posterior 
PDF  P(m  |  d))  assesses  the  probability  of  a  model  under  certain  model  parameters  ( m  =  jx0,JP0,  V,  tQ  j )  given  a 
specific  dataset  (  d  =  j  t ; ,  6t  j  arrival  times  t.  and  back  azimuth  0. )  and  is  written  as: 

P{m  |  d )  =  c(d)P(m)P(d  \  m)  (l) 

Where,  P{rri) is  the  prior  PDF,  P(d  \  m)  the  likelihood  PDF  and  c(d )  a  normalization  function.  P(d  \  m)  is  the 
joint  distribution  of  the  errors  in  back  azimuth  0.(^.  |  m)and  origin  time  O  .{Q.  \  m)  between  the  data  and  the 
forward  model  evaluated  in  each  point  of  the  parameter  space  (Modrak  et  al.,  2010): 


P(d  |  m)=Y\®i{Oi  |  m)d)i(Oi  \  m) 

i= 1 


(2) 
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where,  fl  is  the  number  of  stations  in  the  network,  0  (0.  I  m)  =  .  exp 

■  ■  ppp. 


i  <  r,  V 


and 


the  error  distribution  of  back  azimuths  and  arrival  times,  G  e  and  G ^ 


the  corresponding  variances,  y.  =  9i  —  9p  8i  —  /  —  tp  Oi  and  t(  the  back  azimuth  and  arrival  time  at  station  i  of 
a  candidate  source  of  the  parameter  space.  P(jri)  accounts  for  potential  prior  information  that  can  be  used  to 
constrain  the  parameters.  The  prior  distribution  can  be  written  as: 


P{m)  =  p(x0  ,v0  )p(t0  )p(v)  (3) 


Marginalization 

The  posterior  PDF  ((1)  can  be  marginalized  over  any  dimension  of  the  parameter  space.  BISL  marginalizes  over  the 
velocity  ( V ),  and  time  or  origin  ( t() )  and  obtains  a  new  posterior  PDF  over  the  probable  event  location  jx0 ,  y0  j  : 

^({Wo}  I  d)  =  J  c(d)P(m)P(d  |  m) 

=  P(x0,y0)  {  p(t0)p(v)P[d  I  {x0,j0,^,v})  (4) 

=  c'(d)p(x0,y0)p(d  I  {x0,j0}) 

where,  cr(d )  is  the  new  normalization  function.  The  marginalization  process,  integration  over  specific  dimensions, 
allows  the  assimilation  of  prior  information  in  the  form  of  the  prior  PDF  ((3).  Earlier  implementations  of  BISL 
(Arrowsmith  and  Whitaker,  2008)  use  very  basic  priors  that  just  constrain  upper  and  lower  limits  of  the  sound 
propagation  (0.35  to  0.22  km/s,  respectively)  and  assign  flat  probabilities  to  the  velocities  within  this  range.  This 
velocity  range  is  selected  to  assess  group  velocity  of  rays  refracting  in  the  thermosphere  (0.22-0.24  km/s), 
stratosphere  (0.28-0.31  km/s),  and  troposphere  (0.30-0.34  km/s)  (Ceplecha  et  al.,  1998).  The  approach  taken  here  is 
to  create  prior  PDFs  for  group  velocity  that  assign  different  probabilities  to  sub-regions  of  the  0.22-0.35  km/s 
interval.  By  performing  numerical  simulations  of  infrasound  propagation  using  a  comprehensive  set  of  atmospheric 
scenarios,  we  create  propagation  catalogs  that  are  used  to  evaluate  the  probability  of  occurrence  of  group  velocity 
for  a  specific  range  (  ra ),  azimuth  ( (j)a ),  and  time  (ta)  such  that: 

p(v)  =  K(v,ra,<f>a,ta)  (5) 

where,  K  is  a  function  that  retrieves  and  normalizes  information  from  the  propagation  catalogs  such  that 

K{v,  Va ,  (j)a ,  t  )dv  =  1 .  Assuming  that  there  is  not  prior  information  for  the  location  and  time  of  origin  of  the 

J  V 

event,  we  can  write p(xQ,y0)p(t0)  =  1  .  Using  (5,  (3  can  be  written  as: 

P(m)  =  K(v,ra,<f>a,ta)- 


Networks  (groups  of  sensor  arrays)  can  have  complete  or  partial  azimuthal  coverage  while  recording  an  event. 
Networks  with  arrays  located  at  similar  distances  and  azimuths  may  see  similar  atmospheric  conditions,  and 
therefore  a  unique  prior  PDF  can  be  assigned  to  them.  On  the  other  hand,  networks  with  good  azimuthal  coverage 
may  see  very  different  atmospheric  conditions  depending  on  the  highly  anisotropic  wind  field.  In  the  latter  case,  a 


723 


2012  Monitoring  Research  Review:  Ground-Based  Nuclear  Explosion  Monitoring  Technologies 


unique  prior  PDF  based  on  a  single  azimuth,  range,  and  time  is  not  adequate  to  provide  information  to  the  entire 
network.  Below  we  discuss  future  implementations  of  BISL  that  will  properly  account  for  such  scenarios.  For  now, 
a  more  appropriate  prior  PDF  can  be  formed  by  averaging  the  prior  PDFs  for  each  array  of  the  network  into  a 
general  prior  PDF : 

n 

P(m)  =  kYJKi(v,ra,</>a,ta )  (6) 

i= 1 


where,  H  is  the  number  of  arrays  in  the  system,  and 


dv 


a  new  normalization  coefficient. 


Azimuth  Dependent  Celerity-Range  Histograms  and  Prior  PDF  Generation. 

We  generated  two  propagation  catalogs  using  model-based  regional  and  global  atmospheric  specifications. 
Infrasound  can  propagate  at  local  (few  10s  km),  regional  (100s  km),  and/or  global  distances  (1000s  km).  At  very 
local  (few  kilometers)  and  local  distances,  topography  and  source  geometry  can  play  a  significant  role  in 
propagation.  On  the  other  hand,  propagation  at  regional  and  global  distances  is  mainly  driven  by  the  atmospheric 
structure.  These  specifications  are  generated  by  the  NRL-G2S  atmospheric  model  (Drob  et  al.,  2003)  and  include  air 
temperature,  and  zonal  and  meridional  wind  from  the  ground  to  140  km.  The  regional  dataset  comprises  atmospheric 
specifications  every  6  hours  for  eleven  years  for  a  location  centered  at  the  Utah.  The  global  dataset  comprises 
atmospheric  specifications  taken  within  a  one-year  period  for  locations  randomly  distributed  around  the  entire 
surface  of  the  Earth.  Ray  tracing  uses  the  tau-p  method  (Garces  et  al.,  1998;  Drob  et  al.,  2003)  for  rays  launched  at 
azimuths  covering  360°  in  steps  of  30°,  and  elevation  angles  between  1°  and  50°  at  2°  interval  and  extending  out  to 
2000km  and  12000  km  horizontal  distance  for  regional  and  global  catalogs.  As  a  result  of  these  simulations,  we 
extract  celerity,  maximum  elevation,  and  range  for  multiple  bounces  for  rays  that  are  refracted  to  the  ground.  To 
display  these  results,  we  introduce  the  Celerity  Range  Histogram  (CRH),  which  displays  the  ratio  between  the 
number  of  rays  that  reach  the  ground  and  the  number  of  profiles  at  specific  ranges  and  celerities.  The  CRHs  use: 

1)  an  evenly  discretized  range,  from  0  to  2000  km  and  0  to  10000  km  for  regional  and  global  catalogs  respectively, 
with  50  km  interval,  and  2)  evenly  spaced  celerities  from  0.2  to  0.35  km/s  with  0.01  km/s  interval. 

As  the  propagation  of  acoustic  energy  is  controlled  by  the  highly  anisotropic  vertical  effective  velocity  profiles,  we 
developed  azimuth  dependent  CRHs  by  grouping  ray-path  features  in  four  groups  (quadrants)  depending  on  launch 
azimuth.  Quadrants  covering  315-45  and  135-225  degrees  are  used  to  capture  the  influence  of  meridional  winds  in 
propagation.  The  45-135  and  225-315  quadrants  are  expected  to  capture  the  influence  of  zonal  winds.  In  order  to 
capture  the  variability  in  time  of  the  propagation,  mainly  related  to  changes  in  the  stratospheric  winds  (Whitaker  and 
Mutschlecner,  2008),  different  CRHs  were  computed  for  different  parts  of  the  year.  Figure  1  shows  two  regional 
azimuth  dependent  CRHs  for  the  summer  and  winter  of  2005.  The  two  main  components  in  the  histograms 
(0.29-0.32  km/s  and  0.22-0.26  km/s)  correspond  to  the  expected  thermospheric  and  stratospheric  returns.  Seasonal 
variability  (summer-winter)  is  noted  specially  for  stratospheric  returns  at  azimuths  influenced  by  zonal  winds. 
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Figure  1.  Azimuth  dependent  Celerity  Range  Histograms  for  Winter  2005-2006  (panel  a)  and  Summer  2005 
(panel  b).  CRHs  display  the  ratio  (color  axis)  between  the  number  of  rays  that  return  to  ground  and 
the  associated  number  of  profiles  (used  in  the  propagation  modeling)  at  specific  celerities  (y-axis) 
and  ranges  (x-axis). 
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Figure  2  shows  azimuth  dependent  CRHs  for  the  global  catalog.  It  is  interesting  to  notice  that  for  distances  above 
1000  km  the  histograms  are  very  homogenous  and  seasonal  variation  is  weaker  than  the  effect  for  regional  distances. 
This  effect  can  be  related  to  the  random  distribution  of  the  locations  used  to  get  the  specifications,  as  the  seasonal 
wind  patterns  in  northern  and  southern  hemispheres  have  opposite  behaviors. 
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Figure  2.  Global  azimuth  dependent  CRH. 
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Using  CRH-based  Prior  PDF  to  Enhance  BISL. 

Every  point  of  the  jx0,  y0  j  parameter  sub-space  has  CRH-based  prior  PDFs  for  each  element  of  the  network. 

Station-specific  prior  PDFs  are  averaged  into  one  celerity  distribution  that  is  used  to  marginalize  the  celerity  axis. 
We  tested  the  performance  of  BISF  with  both  uniform  and  CRH-based  prior  PDFs  using  a  surface  explosion.  This 
event  was  a  100-ton  ANFO  explosion  conducted  at  the  Sayarim  Military  Range,  Negev  desert  in  Israel  on  January 
26,  2011  (Gitterman  et  al.,  2011).  The  event  was  recorded  for  three  stations  of  the  International  Monitoring  System 
(IMS)  in  I31KZ,  I46RU,  and  I34MN,  with  distances  of  up  to  6300  km.  Figure  3a  shows  the  localization  of  the 
detecting  IMS  stations  and  Sayarim,  this  figure  also  shows  the  results  of  using  flat  and  CRH-based  prior  PDFs  for 
localization  using  BISF.  These  results  show  a  significant  enhancement  in  the  localization  precision  (decrease  in  the 
95%  credibility  area).  However,  both  results  display  an  offset  between  estimated  and  actual  event  location  (red  star) 
that  can  be  related  to  the  influence  of  wind  in  long-range  propagation.  In  the  current  implementation  of  the  BISF 
algorithm  back  (  0)  and  launch  ( (fi )  azimuths  correspond  geometrically,  such  that  6  =  (/)+  7T .  However,  wind  can 
introduce  propagation  offsets  perpendicular  to  the  launch  angle  and  change  the  expected  back  azimuth.  We  are 
exploring  optimal  techniques  to  correct  for  this  effect. 


IMS  s  unions  (blue  square).  95%  cnxNllbiillty  flal  prtor(green) 
„  95%  credrtibflily  dynamic  prltHmsgerHa) 
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SM  150  WO  350  300  350  «G  ttO 

Speed  (m/9) 


Figure  3.  Using  uniform  and  CRH-based  prior  PDFs  for  localization  using  BISL.  a)  Location  of  the  three 
stations  of  the  IMS  (blue  squares)  that  recorded  and  identified  the  Sayarim  event  (red  star),  95% 
credibility  areas  using  flat  (green)  and  CRH-based  (magenta)  priors,  b)  NRL-G2S  based 
atmospheric  profiles  at  the  location  of  the  explosion. 

Analytical  Expression  for  Model-Based  CRH 

We  are  exploring  the  generation  of  analytical  CRHs  based  on  the  decomposition  of  model-based  CRHs.  We  form 
vertical  cross-sections  of  an  analytical  CRH  (celerity  distributions  for  specific  ranges,  azimuths,  and  times)  by  linear 
combination  of  three  probability  density  functions  (PDFs)  each  characterizing  thermospheric  (fH), 

stratospheric  ( fs )  and  tropospheric  ( fT )  celerity  distributions: 


VtIt  (Mt  ’  ® T  )  Vsfs^S’*  ® s)  'HhIh  ® H  )’  ^ 

where,  /ux  and  CJX  are  mean  and  standard  deviation  of  an  interval  X,  and  ?jx  a  weight  coefficient  for  the  same 
function  in  (0,1)  and  the  sum  of  weight  functions  equal  to  one.  (7  is  commonly  known  as  a  mixed  PDF.  This 
decomposition  of  elements  is  intended  to  provide  a  physical  view  of  the  CRHs,  which  can  be  interpreted  as  the 
superposition  of  the  three  main  atmospheric  ducts.  The  parameters  and  coefficients  required  for  this  decomposition 
are  estimated  using  model-based  CRHs.  First,  the  celerity  is  regrouped  in  three  intervals:  1)  0.22-0.27  km/s, 
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2)  0.28-0.31  km/s,  and  3)  0.30-0.34  km/s.  Second,  the  hit/profile  ratios  for  these  new  intervals  are  added,  and  mean 
and  standard  deviation  are  estimated.  These  nine  parameters  (three  for  each  interval)  are  passed  to  the  (7  to  generate 
'P  which  is  the  analytic  version  of  K .  The  function/*  can  have  various  distribution  shapes.  We  used  the  normal 

distribution  in  this  study;  however,  further  analysis  should  determine  the  most  appropriate  shape  distribution  for  f . 
Figure  4  shows  the  analytic  version  of  the  CRH  of  Figure  lb.  The  intensity  and  width  of  the  thermospheric  and 
stratospheric  returns  are  well  recovered  in  this  analytic  CRH  compared  to  the  model-based  one. 
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Figure  4.  Analytic  version  of  the  CRH.  This  figure  is  the  analytic  version  of  the  summer  version  of  the  CRH 
showed  in  Figure  1. 


CONCLUSIONS  AND  RECOMMENDATIONS 


Implementing  azimuth  and  range  dependent  prior  information  for  group  velocity  in  the  BISL  algorithm  can 
significantly  enhance  the  precision  of  event  localization.  The  azimuth  dependent  CRHs  prove  to  be  an  effective  tool 
to  display  and  summarize  the  results  of  ray  propagation  modeling.  Extracting  and  normalizing  cross-sections  of  the 
CRH  can  easily  generate  prior  PDFs  for  specific  ranges.  By  using  CRH-based  prior  PDFs  we  were  able  to  locate  an 
event  and  determine  a  95%  confidence  area  that  is  significantly  smaller  than  the  one  using  flat  priors.  Our  ongoing 
research  is  focused  now  on  studying  further  enhancements  in  the  application  and  generation  of  prior  information. 
Several  mathematical  representations  of  the  weighting  proportions  Tjx  in  Equation  7  are  under  development. 

The  implementation  of  the  BISL  algorithm  described  here  assumes  a  uniform  CRH-based  prior  PDF  for  the  entire 
network  (averaged  version  of  the  array-specific  prior  PDFs).  However,  even  at  regional  distances  each  station  of  the 
network  has  different  pictures  of  the  atmosphere  between  source  and  receivers.  This  BISL  implementation  cannot 
accommodate  for  station  specific  prior  information  as  it  sees  the  network  as  a  whole  entity  in  the  4D  parameter 
space.  Further  work  would  require  the  manipulation  of  the  BISL  mathematical  framework  to  include  station-specific 
priors.  This  work  may  require  probing  more  complex  parameter  space  that  can  include  expanding  the  velocity  from 
a  parameter  to  a  vector  of  parameters,  with  one  velocity  parameter  for  each  array.  This  expansion  would  allow 
applying  array- specific  prior  PDFs,  however,  this  will  increase  the  dimensions  of  the  parameter  space  and  the 
corresponding  marginalization  can  slow  the  application  of  the  algorithm.  Our  on-going  research  also  includes  code 
development  for  Monte  Carlo  integration  of  the  BISL  posterior  kernel  to  efficiently  probe  these  more  complex 
spaces  and  to  enable  real-time  or  near  real-time  applications. 


727 


2012  Monitoring  Research  Review:  Ground-Based  Nuclear  Explosion  Monitoring  Technologies 


ACKNOWLEDGEMENTS 


We  would  like  to  acknowledge  David  Green  (AWE  Blacknest)  and  Alexis  Le  Pichon  (CEA)  for  the  helpful 

discussions  that  motivated  the  approach  taken  in  this  document. 

REFERENCES 

Arrowsmith,  S.  J.  and  R.  Whitaker  (2008).  Inframonitor:  A  tool  for  regional  infrasound  monitoring,  in  Proceedings 
of  the  2008  Monitoring  Research  Review  Proceedings:  Ground-Based  Nuclear  Explosion  Monitoring 
Technologies ,  LA-UR-08-05261,  Vol  2,  pp.  837-843. 

Ceplecha,  Z.,  J.  Borovicka,  W.  Elford,  D.  ReVelle,  R.  Hawkes,  V.  Porubcan,  and  M.  Simek  (1998).  Meteor 
phenomena  and  bodies.  Space  Science  Reviews  84(3),  327-471. 

Drob,  D.  P.,  J.  M.  Picone,  and  M.  Garces  (2003).  Global  morphology  of  infrasound  propagation.  J.  Geophys.  Res. 
108(D21),  4680. 

Garces,  M.  A.,  R.  A.  Hansen,  and  K.  G.  Lindquist  (1998).  Traveltimes  for  infrasonic  waves  propagating  in  a 
stratified  atmosphere,  Geophys.  J.  Int.  135(1),  255-263. 

Gitterman,  Y.,  J.  Given,  J.  Coyne,  R.  Waxier,  J.  Bonner,  L.  Zerbo  and  R.  Hofstetter  (2011).  Large-scale  controlled 
surface  explosions  at  sayarim,  israel,  at  different  weather  patterns,  for  infrasound  calibration  of  the 
international  monitoring  system,  in  Proceedings  of  the  2011  Monitoring  Research  Review:  Ground-Based 
Nuclear  Explosion  Monitoring  Technologies ,  LA-UR-1 1-04823,  Vol.  2,  pp.  766-777 . 

Modrak,  R.  T.,  S.  J.  Arrowsmith,  and  D.  N.  Anderson  (2010).  A  bayesian  framework  for  infrasound  location. 
Geophys.  J.  Int.  181(1),  399^-05. 

Whitaker,  R.  W.  and  J.  P.  Mutschlecner  (2008).  A  comparison  of  infrasound  signals  refracted  from  stratospheric  and 
thermospheric  altitudes.  J.  Geophys.  Res.  113(D8),  D08117. 


728 


